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Interactions among Arbitrary Materials in Arbitrary Geometries 

M. T. Homer Reid^, Jacob White^'^, and Steven G. Johnson^'-^ 

^Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 

^Department of Electrical Engineering and Computer Science, 

Massachusetts Institute of Technology, Cambridge, MA 02139, USA 

"^ Department Of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 

(Dated: February 24, 2013) 

This paper presents a new method for the efficient numerical computation of Casimir interactions 
between objects of arbitrary geometries, composed of materials with arbitrary frequency-dependent 
electrical properties. Our method formulates the Casimir effect as an interaction between effective 
electric and magnetic current distributions on the surfaces of material bodies, and obtains Casimir 
energies, forces, and torques from the spectral properties of a matrix that quantifies the interactions 
of these surface currents. The method can be formulated and understood in two distinct ways: (1) 
as a consequence of the familiar stress-tensor approach to Casimir physics, or, alternatively, (2) as a 
particular case of the path-integral approach to Casimir physics, and we present both formulations in 
full detail. In addition to providing an algorithm for computing Casimir interactions in geometries 
that could not be efficiently handled by any other method, the framework proposed here thus 
achieves an explicit unification of two seemingly disparate approaches to computational Casimir 
physics. 
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I. INTRODUCTION 

This paper presents a new method for the efficient nu- 
merical computation of Casimir interactions between ob- 
jects of arbitrary geometries, composed of materials with 
arbitrary frequency-dependent electrical properties. Our 
method formulates the Casimir effect as an interaction 
between effective electric and magnetic currents on the 
surfaces of material bodies, and obtains Casimir ener- 
gies, forces, and torques from the spectral properties of a 
matrix that quantifies the interactions of these surface 
currents. Our final formulas for Casimir quantities — 
equations M below — may be derived in two distinct 
ways: (a) by integrating the Maxwell stress tensor over a 
closed bounding surface, as is commonly done in purely 
numerical approaches to Casimir computation [T], but 
with the distinction that here we evaluate the surface 
integral analytically; or, alternatively, (b) by evaluat- 
ing a path-integral expression for the Casimir energy, 
as is commonly done in quasi-analytical approaches to 
Casimir physics [2] , but with the distinction that here we 
are not restricted to the use of geometry-specific special 
functions. In this paper, we present these two distinct 
derivations of our master formulas M and compare our 
new approach to existing computational Casimir meth- 
ods. A free, open-source software package implementing 
our method is available [3]; the technical details of this 
and other numerical implementations of our method will 
be discussed elsewhere. 

Results obtained using our new technique have ap- 
peared in print before [IHH], and Refs. [H [S] briefly 
sketched the path-integral derivation of our method, but 
omitted many details. The purposes of the present pa- 
per are to furnish a complete presentation of the path- 



integral derivation and to present the alternative stress- 
tensor derivation, which has not appeared in print before. 
By arriving at identical formulas — our master formulas, 
equation (fTl) — from the two seemingly disparate starting 
points of path integrals and stress tensors, we explicitly 
demonstrate the equivalence of these two formulations of 
Casimir physics. 

In particular, our demonstration of this equivalence 
furnishes an alternative demonstration that the Maxwell 
stress tensor in dispersive media — questionable under or- 
dinary circumstances — is in fact valid in the thermody- 
namic context, as has been argued on other grounds by 
Pitaevskii ^, and by Philbin [TOIIII] in the context of the 
canonical quantization of macroscopic electromagnetism. 
An algebraic equivalence similar to ours, but relating a 
path-integral expression to the energy density instead of 
the Maxwell stress tensor, was demonstrated in Ref. |I2) . 
which used this equivalence to explain why the dispersive 
energy density (which is valid in ordinary electrodynam- 
ics only for negligible dissipation |13j ) is the appropriate 
quantity to consider in the context of thermal and quan- 
tum fluctuations. Our work does for the stress tensor 
what Ref. [12 does for the energy density. (An alter- 
native approach to relating the stress-tensor picture to 
the energy viewpoint was suggested in Ref. [H], but de- 
tails were omitted; also, the method was restricted to ge- 
ometries that admit a separating plane between objects, 
whereas the method of this paper has no such restriction 
and is applicable even to geometries containing objects 
with interpenetrating features.) 

Although Casimir physics has been with us for some 
seven decades [15], the past fifteen years have witnessed 
a renaissance of interest in the field, driven by labora- 
tory observations of Casimir phenomena in an increas- 



ingly complex variety of geometric and material config- 
urations [T6H2O] . Whereas the theoretical methods used 
in the original Casimir prediction [TS] were restricted to 
the case of simple geometries and idealized materials, re- 
cent experimental progress has spurred the development 
of theoretical techniques for predicting Casimir forces 
among bodies of arbitrary shapes and material proper- 
ties. Such general-purpose Casimir methods have typi- 
cally pursued one of two general strategies. 

A first approach pi f2TH26] seeks to exploit geometri- 
cal symmetries by approximating Casimir quantities as 
expansions in special functions (solutions of the scalar or 
vector Helmholtz equation in various coordinate systems) 
that encode global geometric information in a concise 
way. (Techniques of this sort are often known as spectral 
methods [H].) These methods have the virtue of yield- 
ing compact expressions relating Casimir energies, forces, 
and torques to linear-algebra operations (matrix inverse, 
determinant, and trace) on matrices describing the in- 
teractions of the global basis functions. As is commonly 
true for spectral methods, the expressions are rapidly 
convergent (in the sense of obtaining accurate numerical 
results with low-dimensional truncations of the matrices) 
for highly symmmetric geometries, but less well-suited 
to asymmetric configurations, where the very geomet- 
ric specificity encoded in the closed-form Helmholtz so- 
lutions becomes more curse than blessing and requires 
the special-function expansions to be carried out to high 
orders for even moderate numerical precision. 

An alternative approach is a numerical implementation 
of the stress-tensor formulation of Casimir physics pio- 
neered by Lifshitz et al. [Jf, 29 . Here the Casimir force 
on a body is obtained by integrating the Maxwell stress 
tensor — suitably averaged over thermal and quantum- 
mechanical fluctations — over a fictitious bounding sur- 
face surrounding the body; in modern numerical ap- 
proaches [U I3OH32] the integral is evaluated by numer- 
ical cubature (that is, as a weighted sum of integrand 
samples) with values of the stress tensor at each cuba- 
ture point computed by solving numerical electromag- 
netic scattering problems. As compared to the special- 
function approach, this technique has the virtue of great 
generality, as it allows one to take advantage of the wide 
range of existing numerical techniques for solving scat- 
tering problems involving arbitrarily complex geometries 
and materials. The drawback is that the spatial inte- 
gral over the bounding surface adds a layer of conceptual 
and computational complexity that is absent from the 
special-function approach. 

In this paper we show how the best features of 
these two approaches may be combined to yield a new 
technique for Casimir computations. Our fluctuating- 
surf ace- current (FSC) approach expresses Casimir ener- 
gies, forces, and torques among bodies of arbitrary ge- 
ometries and material properties in terms of interactions 
among effective electric and magnetic currents flowing 
on the object surfaces. The method borrows techniques 
from surface-integral-equation formulations of electro- 



magnetic scattering [33] to represent objects entirely in 
terms of their surfaces — thus retaining the full flexibil- 
ity of the numerical stress-tensor method in handling ar- 
bitrarily complex asymmetric geometries — but bypasses 
the unwieldy numerical cubatures of the usual stress- 
tensor approach to obtain Casimir energies, forces, and 
torques directly from linear- algebra operations (matrix 
inverse, determinant and trace) on matrices describing 
the interactions of the surface currents, thus retaining 
the conceptual simplicity and computational ease of the 
usual special-function approach. 

The FSC formulas for the zero-temperature Casimir 
energy, force, and torque are 
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where the precise form of the matrix M(^) is given in 
SectionlTTj Readers familiar with scattering-matrix meth- 
ods for Casimir computations [H [23]~26J will note the 
striking similarity of our equation (flk) to the Casimir 
energy formulas reported in those works (such as equa- 
tion 5.13 of Ref. PI); in both cases, the Casimir energy 
is obtained by integrating over the imaginary frequency 
axis, with an integrand expressed as a ratio of matrix 
determinants. The difference lies in the meaning of the 
matrices in the two cases; whereas the matrix in typi- 
cal scattering-matrix Casimir methods descibes the in- 
teractions of incoming and outgoing wave solutions of 
Maxwell's equations, the matrix in our equations (fl]) de- 
scribes the interactions of surface currents flowing on the 
boundaries of the interacting objects in a Casimir geom- 
etry. This distinction has important ramifications for the 
convenience and generality of our method. 

In traditional scattering-matrix Casimir methods, the 
matrix that enters equations like (IT]) describes interac- 
tions among the elements of a basis of known solutions 
of Maxwell's equations propagating to and from the in- 
teracting bodies. Such treatments afford a highly effi- 
cient description of scattering in the handful of geome- 
tries for which analytical solutions are available — such 
as incoming and outgoing spherical waves for spherical 
scatterers, left- and right-traveling plane waves for pla- 
nar geometries, cylindrical wave for cylinders, etc. — but 
may be particularly inefficient for describing more gen- 
eral objects, as, for example, if one attempts to describe 
scattering from a cube using a basis of spherical waves. 
Moreover, practical implementations of these methods re- 
quire significant retooling to accommodate new shapes of 
objects; if, for example, having formulated the method 
for spheres, one wishes instead to treat spheroids, one 
must recompute Maxwell solutions in a new coordinate 
system and reformulate the matrices in equations like (fl]) 
to describe the interactions of these new solutions. 



In contrast, the matrix in our equations (IT]) describes 
the interactions of surface currents flowing on the surfaces 
of the interacting objects in a Casimir geometry, as dis- 
cussed in detail in Section |Tl] A crucial advantage of this 
description is that the basis we use to represent surface 
currents is arbitrary; the basis functions are not required 
to solve the wave equation or any other equation, and 
the choice of basis is thus liberated from the underlying 
physics of the problem — we are free to choose a basis that 
efficiently represents any given geometry. One convenient 
choice — though by no means the only possibility — is a 
basis of localized functions conforming to a nonuniform 
surface-mesh discretization (Figure II A ) , where the mesh 
may be automatically generated for arbitrarily complex 
geometries [M] . A particular advantage of this type of ba- 
sis is that, once we have implemented our method using 
basis functions of this type, we can apply it to arbitrary 
geometries with almost no additional effort; in particu- 
lar, having applied the method to spheres, it is essentially 
effortless to apply it to cubes (Section VI I. 

The objective of this paper is to provide two sepa- 
rate derivations of the master formulas (IT]) , one based on 
the stress-tensor formalism and making no reference to 
path integration, and a second based on path integrals 
and making no reference to stress tensors. These deriva- 
tions contain a number of theoretical innovations beyond 
the practical usefulness of the method itself; in particu- 
lar, in the stress-tensor derivation we state and prove a 
new integral identity involving the homogeneous dyadic 
Green's functions of Maxwell's equations (Appendix [b]) , 
while in the path-integral derivation we introduce a new 
surface- current representation of the Lagrange multipli- 
ers that constrain functional integrations over the elec- 
tromagnetic field, which we expect to be a tool of general 
use in quantum field theory. 

The Casimir method described in this paper is closely 
related to the surface-integral-equation (SIE) formulation 
of classical electromagnetic theory. Although well-known 
in the engineering literature [33j . this technique has not 
been extensively discussed in the Casimir literature, and 
for this reason we begin in Section ITT] with a brief review 
of SIE theory. Although the majority of this section is a 
review of standard material, the explicit expressions for 
dyadic Green's functions that we derive in |IIC| have not, 
to our knowledge, appeared in print before. In Sections 
|III| and |IV[ which constitute the centerpiece of the pa- 
per, we present two separate derivations of the master 
FSC formulas (T]); one derivation starts from the stress- 
tensor approach to Casimir physics (Section III), while 



an independent derivation starts from a path-integral ex- 
pression for the Casimir energy (Section IV). In Section 



\V\ we note an important practical simplification that fol- 
lows from the structure of the matrices in equations (JTl). 
In Section |Vl] we validate our method by using it to re- 
produce known results, then illustrate its generality by 
applying it to geometries that would be difficult to ad- 
dress using existing Casimir methods. (Further examples 
of the utility of our method may be found in Refs. jlHH]-) 



Our conclusions are presented in Section [VII[ and a num- 
ber of technical details are relegated to the Appendices. 
Technical details of practical numerical implementations, 
as well as additional computational applications, will be 
discussed elsewhere. 



II. A REVIEW OF THE 

SURFACE-INTEGRAL-EQUATION 

FORMULATION OF CLASSICAL 

ELECTROMAGNETISM 

Computational Casimir physics is intimately related 
to the theory of classical electromagnetic scattering, and 
many practical methods for predicting Casimir interac- 
tions are based on well-known techniques for solving scat- 
tering problems. Among the classical scattering methods 
that have been appropriated for Casimir purposes are 
the T-matrix method PH [351 [3S] , the method of reflec- 
tion coefficients |53j, and the numerical finite-difference 
method [30l [31]. The method discussed in this paper 
derives from yet another well-known approach to scat- 
tering problems, namely, the method of surface integral 
equations (SIEs). (SIE techniques were first used for 
Casimir studies in Ref. 4|, while Refs. |32l[37| presented 
an SIE-based implementation of the numerical stress- 
tensor method.) As background for the remainder of this 
paper, in this section we review the well-known SIE pro- 
cedure. 

Surface-integral techniques have a long history in elec- 
tromagnetic theory [3^, dating back to the equivalence 
principles of Love [35] and Schelkunoff [ID] and the 
Stratton-Chu equations [41] from the first half of the 20th 
century. Numerical implementation of SIEs — known as 
the "boundary-element method" (BEM) or the "method 
of moments" — emerged in the 1970s as an alternative 
to other computational procedures such as the finite- 
difference method (FDM) and the finite-element method 
(FEM) m US]. Whereas the FDM and the FEM pro- 
ceed by numerically solving a spatially local form of 
Maxwell's equations, and thus allow treatment of mate- 
rials with essentially arbitrary spatial variation of the di- 
electric permittivity and magnetic permeability, the SIE 
approach takes advantage of the analytically known solu- 
tions of Maxwell's equations in homogeneous media, and 
thus in practice is most readily applicable to piecewise- 
homogeneous material configurations. For this reason, 
while the FDM and FEM have the advantage of be- 
ing able to treat a wider class of materials, the SIE 
method exhibits significant practical advantages for the 
piecewise-homogeneous geometries typically encountered 
in Casimir studies. 

To fix ideas and notation for the remainder of this 
paper, we here review the SIE formulation of electro- 



magnetic scattering problems, beginning in Section II A 
with the simplest case of perfectly electrically condu cting 
(PEC) scatterers, and then generalizing in Section II B 
to the case of arbitrary materials. 




FIG. 1. Schematic depiction of a scattering geometry in the surface-integral-equation picture. A collection of arbitrarily-shaped 
homogeneous bodies, with frequency-dependent relative electrical properties {e'^,^'"}, is embedded in a homogeneous medium 
with electrical properties {e'', jj,'^}. Incident radiation, characterized by electric and magnetic fields E™'^, H'"'^, impinges on the 
objects to induce surface currents; for perfectly conducting objects we have only electric surface currents (K), while for general 
objects we have equivalent electric and magnetic (N) surface currents. The goal of surface-integral-equation methods is to solve 
for the surface-current distributions in terms of the incident fields, after which we can compute the scattered fields anywhere 
in space from the surface currents. [The dotted line indicates a fictitious bounding contour C surrounding one of the objects 
over which we integrate the Maxwell stress tensor to compute the Casimir force on that object (Section III )] 



The material of these two subsections is well-known 
and entirely standard within the computational electro- 
magnetism literature, and is reviewed here only for com- 
pleteness. However, in Section II C we extend the SIE 



formalism one step beyond what is usually done to write 
explicit expressions [equations (20) and (23l] for scatter- 



ing dyadic Green's functions in terms of the SIE matrices 
and the homogeneous dyadic Green's functions (DGFs). 
Although these expressions are straightforward conse- 
quences of the standard BEM procedure outlined in Sec- 
tions IIA-B, to our knowledge they are appearing here 
for the first time. 

Throughout this section we will refer to the scattering 
situation depicted schematically in Figure [T] in which 
a collection of homogeneous scatterers (with frequency- 
dependent relative electrical properties {e'",/i''}) is em- 
bedded in a homogeneous medium (electrical properties 
{e'^,/i'^}) and irradiated by incident radiation character- 
ized by an incident electric field E'"'^. 



A. The SIE Method For PEC Bodies 

We first consider the case in which the scattering ob- 
jects in FigurefTlare perfect conductors. An incident field 
impinging on PEC bodies induces a tangential electric 



current distribution K(x) on the body surfaces, which 
gives rise to a scattered field according to 



E""^*(x) = / r™''=(x,x') •K(x')dx'; 



(2) 



here the integral extends over the surfaces of the bod- 
ies and r"^^'*^ is the homogeneous DGF for the exterior 
medium. (Our notation for DGFs is summarized in Ap- 
pendix|Aj throughout this section we work at a single fre- 
quency and suppress frequency arguments to E, F, and 
K.) For a given incident field E'"'^ we can solve for K by 
requiring that the total (incident -I- scattered) field sat- 
isfy the appropriate boundary condition, which for PEC 
bodies is simply that the total tangential E-field vanish 
for all points x on the body surfaces: 



gscat^^^_^gi„c(^^ Xn(x)=0. 



[Here taking the cross product with n(x), the outward- 
pointing surface normal at x, is simply a convenient way 
of extracting the tangential components of a vector.] In- 
serting (pi) yields an integral equation for K(x) : 



'(x,x')-K(x')dx' 



X n(x) = -E'"'=(x) X n(x). 
(3) 
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TABLE I. A glossary of symbols used in this paper. 



Equation ([3]) is known as the "electric-field integral equa- 
tion (EFIE7' 33 . 

Thus far all we have done is to restate the problem in 
an integral-equation form. The next step is to discretize 
this integral equation by introducing a finite set of tan- 
gential vector-valued basis functions {fQ(x)}, defined on 
the surfaces of the bodies, which serve a dual purpose as 
expansion functions for surface currents and test func- 
tions for boundary conditions. As noted in Section IT] an 
advantage of SIE methods is that they place no restric- 
tion on these basis functions; in particular, the {fa} need 
not solve the wave equation or any other equation and 
need not encapsulate any global information about the 
scattering geometry. Of course, if symmetries are present 
then we may wish to choose the {fa} in a way that re- 
flects them — we might choose vector spherical harmonics 
for a spherical scatterer, say, or a Fourier basis for a 
planar scatterer — but nothing in the SIE formulation re- 
quires such a choice, and we are equally free to choose the 
{fo,} to be arbitrary polynomials, piecewise-linear func- 
tions, or any other functions we like. For scatterers of 
complex geometries, a particularly convenient strategy 
is to discretize object surfaces into small flat panels and 
take the {£„} to describe elemental currents sourced and 
sunk at panel vertices [H], as depicted in Figure II A 



The localized basis functions that result from such a pro- 
cedure are known as "boundary elements," and SIE im- 
plementations based on them are commonly known as 
"boundary-element methods" (BEM) or the "method of 
moments." 

Having chosen a set of basis functions, the surface elec- 
tric current distribution is approximated as a finite ex- 



pansion in the {fc}: 

K(x)=^fc„f„(x). 



(4) 



This expansion is then inserted into (|3|, and the inner 
product of that equation is taken with each member in 
the set {fa}, yielding one equation for each of the un- 
known coefficients ka. Collecting these equations yields 
a linear system of the form 



Mk = v, 



(5) 



where k is the vector of ka coefficients, the elements of 
the RHS vector v describe the interactions of the basis 
functions with the incident field. 



/sup ft 



f„(x)-E'"'=(x)dx 



j.inc\ 



/' 



(6a) 
(6b) 



(where supfa is the support of basis function fa), and 
the elements of the M matrix describe the interactions of 
the basis functions with each other through the exterior 
medium: 



Map = / dx / dx' f,(x) • r-^''^(x,x') • f/3(x') 

(7a) 
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fa 



(7b) 



The linear system (l5| may be solved for the surface- 
current expansion coefficients {ka}, after which we can 
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FIG. 2. One possible choice of expansion functions for tangential currents on object surfaces is obtained by discretizing object 
boundaries into unions of small flat panels and introducing localized basis functions describing elemental currents sourced and 
sunk at panel vertices: these are "RWG" basis functions, associated with each edge (pair of triangles) in the mesh [44j . 



compute the components of the scattered field at an arbi- 
trary point in the exterior medium using the discretized 
form of equation (pi): 



^f'^*(x) 



/^' 



r:=f^(x,x')/.,(x')dx' 



sup fo. 



^Y.kJTT'\^) 



fa 



(8b) 



[The second lines of equations (l6|, (It]), and (Is]) define 
some useful shorthand for commonly-encountered inte- 
grals; in (Isb) note that the contracted index and the in- 
tegrated argument to T are suppressed, while the uncon- 
tracted index and non-integrated argument are written 
out.l 



B. The SIE Method For General Bodies 

For PEC bodies, the mathematics of the SIE procedure 
neatly mirrors the physics of the actual situation. Indeed, 
for good conductors at moderate frequencies it really is 
true that the physical induced currents are confined near 
the object surfaces; the surface current distribution K(x) 
for which we solve in an SIE method thus has a direct 
physical interpretation as an induced surface current. 

The situation is more complicated for general (non- 
PEC) objects, for here the physical induced currents are 



no longer confined to the surfaces, but instead extend 
throughout the bulk of the object. The obvious exten- 
sion of the procedure outlined above would be to intro- 
duce a volume discretization and solve a system analo- 
gous to (pi) for the coefficients in an expansion of a volume 
current density J(x). Such a procedure, while retaining 
the intuitive interpretation of the quantity computed as 
a physical current density, would suffer from poor com- 
plexity scaling, as the number of unknowns [and thus 
the dimension of the linear system corresponding to ([5|] 
would scale like the volume, not the surface area, of the 
scattering objects. 

An alternative approach is to abandon the strategy 
of solving for the physical volume sources and to solve 
instead for equivalent surface sources that give rise to 
the same scattered fields. The mathematical machin- 
ery underlying this approach is a vector generalization 
of Green's theorem known as the Stratton-Chu equa- 
tions [41] , which relate the E and H fields in the interior 
of a region to the tangential components of the fields on 
the boundary of that region. More precisely, let dOr be 
the surface of the rth object in our geometry, and for 
points X on dOr define two tangential vector fields ac- 
cording to 



K'='*(x) = n(x) X H(x) 



N'="(x) = E(x) X n(x) (9) 



where n is the outward-pointing normal to dOr at x and 
where E and H are the total fields at that that point. The 



Stratton-Chu equations are then the following surface- 
integral expressions for the fields inside and outside dOr'- 
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(lib) 



In equations (10 11), the r and e superscripts on T 



label the homogeneous DGFs for the medium interior to 
dOr and the exterior medium, respectively; the spatial 
arguments to E, F, K, and N are as in equation (pi), but 
are suppressed here to save space.] 

Note the following differences between expressions (fTol) 
and (11) for the interior and exterior fields: (a) the sur- 
face integrals in the two cases differ in sign, arising from 
the reversal of direction of the surface normal in ^ ; (b) 
the r dyadics in ( 10 ) are those for the homogeneous 
medium interior to Or, while in (11) we instead have 
those for the the exterior medium; (c) in (|10| ) we inte- 
grate over the single surface dOr, while in ( |11[ ) the inte- 
gral is over the union of all object surfaces, U dOr (which 
we may think of as the boundary of the exterior medium, 
UdOr ~ dOe)\ (d) the incident fields contribute to ex- 



pressions ( 11 ) for the exterior fields, but are absent from 
expressions ( 10 ) for the interior fields. 

Although the tangential vector fields defined by ^ 
are simply the components of the E and H fields and 
do not correspond to physical source densities, nonethe- 
less the form of equations (llO 11) suggests interpreting 
'K°^ and N°^ as effective electric and magnetic surface 
current densities, which, if known, would allow computa- 
tion of the fields anywhere in space, just as knowledge of 
the physical surface current K suffices in the PEC case 
to determine uniquely the full scattered field. To em- 
phasize this analogy, we will henceforth drop the "eff" 
designation from K and N. 

As in the PEC case, the K and N distributions are de- 
termined by requiring that the total fields satisfy appro- 
priate boundary conditions. For non-PEC bodies these 
are simply that the tangential components of the total 
fields be continuous across material boundaries; for a 
point X on the surface of a body we have 



E°"'(x)-E'"(x) xn(x)=0 
H°"'(x)-H'"(x)l xn(x) = 0. 



(12a) 
(12b) 



Inserting ( lOpT ) into ( 12 ) leads to integral equations for 
K and N that generalize equation (Is]) for non-PEC bod- 
ies. As in the PEC case, the next step is to discretize 



these integral equations by approximating the electric 
and magnetic surface currents as expansions in a finite 
set of tangential vector- valued basis functions defined on 
the object surfaces, 



K(x) 






fc„f„(x), N(x) = -^n„f„(x), (13) 



and testing the integral equations obtained from ( 12 ) 



with each basis function. (The minus sign in the mag- 
netic surface-current expansion is a useful convention 
that leads to a symmetric linear system [13].) The re- 
sult of this procedure is a linear system of the same gen- 
eral form as (Isl), but now enlarged to exhibit a 2x2 block 
structure: 



jyjEE ypm 
jyjME ]y[MM 





(14) 



In equation (14), the elements of the RHS vector describe 



(15) 



the interactions of the basis functions with the incident 
electric and magnetic fields [compare equation (ro|], 

while the elements of the M matrix describe the basis 
functions interacting with each other both through the 
exterior medium and through the medium interior to one 
of the scattering bodies. For example, the elements of the 
M'^^ block are 






t/3 



(16) 



and similarly for the other blocks. (The T^'^'^ term here 
is present only if basis functions f^ and f^ are defined on 
the surface of the same object Or, while the T^^'"^ term 
is present even for basis functions defined on the surfaces 
of different objects.) 

After solving (14), the scattered fields at an arbitrary 
point X are obtained, in analogy to equation ([8|, from 
the expansions 



i?r(x) = 5]{fc„(rr^(x) 

a 



f„)- 



(X)|fa)} 



(17a) 



(x)|f.)} 



(17b) 

(These are the scattered fields in the exterior region; 
the expressions for fields in the interior of object r are 
similar, but involve the homogeneous DGFs r^''''' for the 
medium interior to object r.) 



Explicit SIE Expressions for Dyadic Green's 
Functions 



The discrctized SIE method reviewed in the previous 
two subsections is typically employed as a numerical tech- 
nique, with the linear systems (p]) and (14) solved using 



methods of computational linear algebra and the scat- 
tered fields (Isl) and (171 evaluated numerically. In this 



paper, in contrast, we will use the SIE formalism in a 
somewhat unusual way, by carrying the analytical de- 
velopment one step further than is commonly done. By 
exploiting the formal solution of equations ^ and (14 1, 
we will obtain useful expressions for scattering dyadic 
Green's functions in terms of the formal inverse of the 
SIE matrix M. These expressions will then be used in 
Section |III| to derive compact FSC expressions relating 
Casimir quantities to linear-algebraic manipulations of 
the M matrix. Although these final expressions will ul- 
timately be evaluated numerically, the analytical expres- 
sions derived in this subsection are an important ingre- 
dient in their derivation by stress-tensor methods. (The 
expressions derived in this subsection are not needed for 
the path-integral derivation of the FSC Casimir formu- 
las.) 



of the scattered electric field, which has the effect of sub- 
stituting r"^ for r'='= in ([2|. The resuh is 



(x,x')--5:(rr^(x 



a/3 



fa VK, 



a;3\I/J 



(20b) 



The General Case 

To obtain explicit expressions for scattering DGFs in 
general geometries, we mimic the procedure followed 
above, but now using the general SIE formalism out- 
lined in Section IIIBI instead of the PEC formalism of 
Sectior II A To compute G^^, we again take the incident 
field to be the field of a unit-strength j-directed point 
electric source at a point x' in the exterior medium, in 
which case the elements of the RHS of equation ( 14 ) are 



The PEC Case 

The scattering dyadic Green's function t^f?^(x, x') is 
the scattered electric field at x due to a point electric 
source at x' (Appendix IaI); here we will need the case 
in which both x and x' lie in the exterior medium. To 
compute this quantity using the SIE technique of Sec- 
tion |II A[ we take the incident field to be the field of a 
unit-strength j-directed point electric current source at 
a point x' in the exterior medium, which is simply 

i?r(x) = r-'^(x,x'). 

Then the elements of the RHS vector in (ISl) are 



^(x') 



(18) 



while the coefficients in the expansion of the scattered 
field may be obtained as the formal solution of (ISl), 



/3 



(19) 



(where W = M ^ is the inverse SIE matrix). Inserting 
(191 and (18) into ([8|, the scattered field at x — which is 



just the scattering DGF we are seeking to compute — is 






(x,x') = -E 



(x)f,)W^„^(f, 



-pEE,e 
J 



(x') 



(20a) 

We will also need the magnetic-magnetic DGF g^^"^^ 
which is the scattered magnetic field due to a point mag- 
netic source. This is obtained in easy analogy to the 
above by (a) taking the incident field to be the field of a 
magnetic point source instead of an electric point source, 
which has the effect of substituting T^" for F^^ in (llS]) ; 



and (b) computing the scattered magnetic field instead 




(21) 



The expansion coefficients that enter into equation (17) 



are given, in analogy to equation (19), by the formal 
solution of ( 14 ) : 













(22) 



Inserting (21) and (22) into (17), and proceeding simi 



larly for the magnetic-magnetic case, then yields the gen- 
eralization of equation ( 20 ) to non-PEC geometries: 



er(x,x') = 



/ (rr^(x)|f„) 
t^U(rr^(x)|f.) 



Wo^P 



iW. 



3 

ME.e 



(23a) 



'(x,x') = 




(f,|r;^-'^(x'))y 

(23b) 



Wo^P 




Equations (20) and (23) are the most important re- 



sults of this section of the paper. The crucial property of 
these expressions is that they present the inhomogeneous 
Green's function in a fuUy-factorized in which factors 
depending on x are separated from those depending on 
x'. In thise sense, equations (20) are similar to Green's- 



function expansions for special geometries commonly en- 
countered in the literature, such as spherical-harmonic 
expansions for spherical geometries or Bessel-function ex- 
pansions for cylindrical geometries [13l[36]; the difference, 
of course, is that ( 20 ) is applicable to arbitrary geome- 



tries, with the geometric information encoded in the W 
matrix and the basis functions {fa}. 



III. STRESS-TENSOR DERIVATION OF THE 
FSC CASIMIR FORMULAS 

The stress-tensor approach to Casimir physics re- 
lates Casimir forces to classical dyadic Green's functions 
(DGFs). This technique was pioneered by Dzyaloshin- 
skii, Lifshitz, and Pitaevskii (DLP) in the 1950s [251 US] 
and has remained an important computational strategy 
ever since [451 I46j : in particular, modern numerical al- 
gorithms for computing Casimir forces between bodies 
of complex geometries have tended to use the stress- 
tensor approach, with values for the relevant DGFs com- 
puted numerically \T\ f3m - l32] . Here, after briefly review- 
ing the formalism relating Casimir forces to DGFs (Sec- 
tion III A I , we will show that the concise SIE expressions 
for the DGFs that we derived in Section III CI afford a 
significant simplification of the usual computational pro- 
cedure. In particular, we show that the surface integral 
of the stress tensor, which in previous work has typically 
been evaluated by numerical cubature, may in fact be 
evaluated analytically for an arbitrary closed surface of 
integration, leaving behind a simple expression relating 
the Casimir force to the trace of a certain matrix. 



A. A Review of Stress-Tensor Casimir Physics 

In the stress-tensor approach, the i-directed Casimir 
force on a body is obtained by integrating the expectation 
value of the Maxwell stress tensor over a closed bounding 
surface surrounding the body: 



J'^ 



d^ 



F^iO, 



(24) 



F,(C)=£(T.,(e;x))n,(x)dx. (25) 

Here the integration surface C may be the surface of the 
body in question or any fictitious closed surface in space 
bounding the body (as in Figure nl), and the expectation 
value is taken with respect to quantum and thermal fluc- 
tuations. The expectation value of T^- is next written 
in terms of the components of the electric and magnetic 
fields, 



Tij) -elEiEj 



2 



l-i{H,Hj 



e{EkEk 



KHkHk 



(26) 



scattering DGFs 
relations read 



at temperature T = 0, the 



i?.(C,x)i?,(e,x'))--;i€e,T(C;x,x') 



i/,(^,x)ij,(^,x')) = -HQlTii:^,^) 



(27a) 
(27b) 



where, as discussed in Appendix \K\ tj'=^^(^,x, x') is the 
scattered portion of the electric field at x due to an elec- 
tric current source at x', all quantities having time depen- 
dence ex e"*"^* ; similarly, Q^" gives the scattered magnetic 
field due to a magnetic current source. (In the original 
work, DLP wrote V x V x g"""^ in place of g""; the 
equivalence of the two quantities has been discussed e.g. 
in Ref. 11 .) 



Inserting (261 and (27) into (25) yields an expression 



for the Casimir force-per-unit-frequency in terms of scat- 
tering DGFs: 



F.iO 



K 



eG^+f^G. 



MM_ 

ij 



2 



^Gkk+I^Gkk 



Y' 



dx. 

(28) 
Equation ( 28 1 is the starting point of many numerical 



Casimir studies, as it reduces the computation of Casimir 
forces to the computational of classical DGFs. In prin- 
ciple, the DGFs in question may be computed using any 
of the myriad available numerical techniques for classical 
scattering problems; to date, numerical Casimir investi- 
gations using both the finite-difference method 1301 131] 
and the discretized SIE method reviewed in Section 
In] [32l [37] have appeared. In these studies, the surface 
integral in (28 1 is evaluated by numerical cubature, with 



the values of the integrand at each cubature point com- 
puted by solving numerical scattering problems. 

Here we proceed in a different direction. Instead of 
taking equation (28) as the jumping-off point for a nu- 



merical investigation, we will continue the analytical de- 
velopment one step further by inserting our explicit SIE 
expressions ( 20 ) and ( 23 ) into equation ( 28 ) and analyz- 



ing the result. As we will see, this step will allow us to 



evaluate the surface integral in (281 analytically, elimi- 
nating the need for numerical cubature and resulting in 
a compact matrix-trace formula for the Casimir force. 

Because the thrust of the argument is easiest to present 
in the simplest case of perfectly electrically conducting 



(PEC) bodies, we begin with that case in Section IIIB 
leaving the treatment of general materials to Section 

[mcj 



B. Stress- Tensor Derivation of FSC Formulas for 
PEC Objects 



[Here it is understood that e — e^e'^ and /i — figfi'^ are 
the (spatially constant) permittivity and permeability of 
the exterior medium at the frequency in question; e'^,/i'^ 
are the dimcnsionless relative quantities.] Finally, the 
fluctuation-dissipation theorem is invoked to relate the 
expectation values of products of field components to 



In Section [TTC] we derived explicit SIE expressions for 
the DGFs that enter into the integrand of (28); for the 



case of PEC scatterers, the relevant expressions are equa- 
tions (20). Our strategy here will be to insert these ex- 



pressions into (28) and analyze the result; to facilitate 



this procedure, it is convenient first to write equations 
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(201 in a slightly different form by (a) expressing the 



four r dyadics in terms of the two G and C dyadics 
(Appendix [A|), and (b) writing out inner products like 



(f |r) explicitly as integrals over the supports of the basis 
function f [compare equations (|6|, Q, and ( [8|]. Then 
the quantities that enter into the integrand of (|28|) are 



eg-(x,x) = -6^ (rr'-^lx) ^o)wU^p rf ^^(x') 



-^lQ^J.''{n''f^Wafi\ I G,fc(r„,x)/Qfe(r„)(ir„ U / /^^(r^) G£j(x, r^) dr;3 



(29a) 



a/3 



Me™(x,x) = -^^ (rr''^(x) io)w^Mp 



a/3 



^(x' 



= +fJ,OfJ.''{K''f ^ Wa/3 < / Cik{ra,x) /afc(ra) dra\ < / //3/?(r,3) Cejix, Tfj) dVfj 

„fl Wsupfa J I Jsupf,3 



(29b) 



a/3 



(Here k*^ = \JWJJ^ ■ ^ is the imaginary wavenumber of the exterior medium, and we have suppressed the dependence 
of the G and C tensors on k*^.) Note that both of these expressions have the same form: a sum over basis functions 
fa and f/3, with a summand involving integrations over the supports of the basis functions. Indeed, equations (29i) 



and (29d) are identical up to the different kernel functions (G or C) that enter into the integrals over basis functions. 
Note also that the variable x, which is the integration variable in the surface integral in (28), appears in (291 only 



through these kernel functions. This implies that, after inserting (29 1 into (28), we will again have a sum of terms 



of this same form — a sum over basis functions, with the summand involving integrals over the basis functions — and, 
moreover, that many of the factors in this summand will be independent of the integration variable x in (28) and 



may thus be pulled outside the surface integral, which will now contain only factors of G and C. The result is 
(Zo - VAio/eo, Z'' - y^W) 



F,{£,) = -^Wap ■ ZoZ^k" dva drplfak{ra)Itki{ra,ri3)f(}i(rp)\ 

"" Q^ Jsupfc J sup f ft 

where, as anticipated, the surface integral is now contained inside the definition of the I kernel: 

I»fc^(ra,r/3) = [k")^ j> {Gjfc(ra,x)Gfj(x,r^) -Cjfc(r„,x)C£j(x,r^) 

--^ G„fc(ra,x)G£™(x,r^) - Gmfc(ra,x)G£m(x, r^) 



(30) 



hjd'x.. 



The fact that W is a symmetric matrix iWaji — Wpa) allows us to rewrite equation (30) to read 

^^(0 = ^ I] W^"/3 • ZoZ'^' I dv^ f drp f /afc(r„)T,fc^(ro, r^)//3,(r/3)| (31) 

where we have defined a symmetrized version of the I kernel: 

Tike{r,r') =lM{r,r')+Iiik{r',r). 

The point of this step is that, as demonstrated in Appendix |B[ the surface integral in the definition of the I kernel 
may be evaluated in closed form, for any topological two-sphere C, with the result 



'0, if r, r' lie both inside or both outside C 

-^-pG/c^(r', r"^) if r, r' lie on opposite sides of C 



(32) 



where, in the second case, r' (r"^) is whichever of r,r' lies in the interior (exterior) of C 
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Armed with the dichotomy (32 1, we can now analyze the quantity in curly brackets in (31 1. Recall that the bounding 
contour C encloses one of the objects in our Casimir geometry; call this object Oi and the remaining objects C'2,3,---- 



Equation (32) then tells us that the curly-bracketed term in (31 1 vanishes except when precisely one of the basis 



functions \Y^, fp} lies on the surface of object Oi. When this condition is satisfied, the integral over basis functions 
in (31) reads 



^sup fa 



supf^ 



r t? 

drp \fak{Ya) T. Gfcf (r^ , T/j) ffn{Yfj) 

L OYfyi 



drat 

But this is nothing but the derivative of the a, (3 element of the SIE matrix pj with respect to a rigid infinitesimal 
displacement of object Oi in the i direction, 



_d_ 



M, 



af3- 



Inserting this into (31), we find that the imaginary- 



frequency-^ contribution to the Casimir force is given 
simply by 



m) 



2n ^^ 



Wa0 



= — Tr [m-1 
27r L 



"A 

d_ 
dr. 



M, 



aP 



M 



(33) 



(where we have recalled the definition W = M^^), and 
inserting this into (24) we obtain the FSC formula for the 
Casimir force, equation (fib). To obtain the FSC formula 
for the Casimir energy, we note that the Casimir force 
on an object is minus the derivative of the energy with 
respect to a rigid displacement of that object; using the 
standard identity 



d_ 
dr,. 



log det M = Tr 



M" 



dri 



M 



and choosing the zero of energy to correspond to the en- 
ergy of the configuration in which all objects are removed 



to infinite separations (for which configuration we denote 
the SIE matrix by Moo), we recover equation (ITk). Fi- 
nally, equation (Tp) follows from taking derivatives with 
respect to a rigid rotation instead of a rigid displacement. 
This completes the stress-tensor derivation of the FSC 
formulas for the case of PEC objects. 



C. Stress- Tensor Derivation of FSC Formulas for 
General Objects 

The derivation of the FSC formulas for general objects 
is now a straightforward generalization of the procedure 
for PEC objects. Again we start with equation (28), and 



again we insert in this equation the factorized expressions 
for scattering DGFs that we derived in Section fllCI the 
difference is that for non-PEC objects we must now use 
the more complicated expressions (23). Mimicing the 



discussion following equations ( 29 1 above now leads to a 



modified version of equation (311 in which the X kernel 
is promoted to a 2 x 2 matrix: 



FM) 







ap ctp 



dvc / drislfakira) 

/sup fa -/sup f^ 



ZoZ^K^I.kiircrp) K''Jipi{ra,rf3) 
~K''J,pe{ra,rp) ■^^T,pi{ra,rp) 

with Tr denoting a 2 x 2 matrix trace and the JJ kernel defined in analogy to I: 

Jiki{r,r') = J^k£{r,r')+Juk{r',r) 




Jike{ra,rf}) = {k'') (t {Gifc(rc,x)C£j(x,r^) + Cifc(rQ,x)G'<?j(x,r^) 

S. 



■H2_ 

2 



G„fc(ra,x)Q™(x, r^) + C'TOfe(rQ,x)Gfo„(x, r^) \ hjdx. 
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Again in analogy to I, the surface integrals in the definition of J may be evaluated in closed form to yield 

'0, if r, r' lie both inside or both outside C 

— — Ckiiv\ r"^) if r, r' lie on opposite sides of C 
or] 



(35) 



derivative of the SIE matrix 



and, armed with (32 1 and (35), it is now easy to identify the integral over basis functions in (34) as nothing but the 



sup fc 



. ZoZ''K''Iikeira,rp) K^'Jipiiva^rp) . 



supf^ 






(36) 



material SIE matrix in equation ( 14 ) 



Inserting ( 36 ) into ( 34 ) now simply reproduces equation ( 33 ) with the M matrix understood to refer to the general 
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IV. PATH-INTEGRAL DERIVATION OF THE interaction matrix that appears in the SIE formulation 

FSC CASIMIR FORMULAS of scattering reviewed in Section In] 



It is remarkable that the path-integral approach to 
Casimir physics, which bears little superficial resem- 
blance to the stress-tensor formalism of the previous sec- 
tion, may nonetheless be used to furnish a separate and 
entirely independent derivation of the same FSC formu- 
las that we derived above using stress-tensor ideas. In 
this section, after first reviewing the well-known formal- 
ism for obtaining Casimir energies from constrained path 
integrals (Section IV A), we present this alternate deriva- 
tion (Section [TVB^ 

The path-integral procedure presented here differs 
from typical path-integral treatments of Casimir phe- 
nomena in at least two ways. First, whereas many au- 
thors write the action for the electromagnetic field in 
terms of the gauge-independent E and B fields [36] , or in 
terms of the four-vector potential A'^ in a way that de- 
pends on a specific choice of gauge (often the "temporal" 
or "Weyl" gauge A° = 07]), here we write the action 
in terms of A^ with a Fadeev-Popov parameter that al- 
lows arbitrary gauge choices; we verify explicitly that the 
Fadeev-Popov parameter is absent from all final physical 
predictions. (This portion of our treatment is similar to 
that of Ref. gS].) 

Second, we introduce a new implementation of the con- 
straint that the path integral extend only over field con- 
figurations satisfying the boundary conditions. Our rep- 
resentation emphasizes the continuity of the tangential 
E and H fields across the surfaces of the objects in a 
Casimir geometry, and the Lagrange multipliers that we 
introduce to enforce the constraints have an attractive 
physical interpretation as surface currents, thus estab- 
lishing a connection to the SIE ideas reviewed above. Af- 
ter integrating out the photon field, we are left with func- 
tional integrals over surface-current distributions, with 
an effective action describing the interactions of these 
currents through the electrical media interior and ex- 
terior to the objects; upon discretization, this action 
turns out to involve precisely the same surface-current- 



A. A Review of Constrained Path-Integral 
Techniques for Casimir Energies 

Path-integral formulations of field-fluctuation prob- 
lems were pioneered by Bordag, Robaschik, and Wiec- 
zorek [35] and by Li and Kardar [321 HO] and have since 
been further developed by a number of authors (see 
[36l l5T] for extensive surveys of recent developments.) 
In this section we review the key steps in this approach. 



Casimir Energies from Constrained Path Integrals 

In the presence of material boundaries, the partition 
function for a quantum field 4> (which may be scalar, vec- 
tor, electromagnetic, or otherwise, but is here assumed 
bosonic) at inverse temperature /? takes the form 



Z0» ^ I 



2?0(t,x) 



(37) 



where the action Sfj is the spacetime integral of the Eu- 
clidean Lagrangian density for the (j) field. 



Sb\4>] = 



hp 



dr I dx C, 



;{'^(^,x)}. 



(38) 



in (37) indicates that this 



and where the notation [• • 
is a constrained path integral, in which the functional 
integration extends only over field configurations satis- 
fying the appropriate boundary conditions at all material 
boundaries. 

If the boundary conditions are time independent and 
the Lagrangian density contains no terms of higher than 
quadratic order in (/) and its derivatives, then it is conve- 
nient to introduce a Fourier series in the Euclidean time 
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variable, 



(T,x)=^0„(x)e-^«"^ Cn 



2Tm 



whereupon the path integral ( 37 1 factorizes into a prod 



uct of contributions from individual frequencies, 



with 



S 



z{P;^n) 



0n(x);^„ 



!?(/.„ (x) 



.-S[0„;?„] 



(39) 



dx/:B{0„(x)e-'«""} 



representing the contribution to the full action ( 38 1 made 



only by those field configurations with Euclidean-time 
dependence ~ e~'^"'^. The free energy is then obtained 
as a sum over Matsubara frequencies, 



/3 Z^{p) 



^ CO 

= --yin 



P 






(40) 



where Zao(Zao) is Z(Z) evaluated with all material ob- 
jects separated by infinite distances [dividing out these 
contributions in ( [40| is a useful convention that amounts 
to a choice of the zero of energy] . In the zero-temperature 
limit, the frequency sum becomes an integral, and the 
zero-temperature Casimir energy is 



h 



d^ln 



2(0 



(Here and below we omit the /3 argument to Z). 

Enforcing Constraints via Functional 5-functions 



(41) 



Equations ( 40|[4T| ) reduce the computation of Casimir 
energies to the evaluation of constrained path integrals 
( 39 1 . In most branches of physics, the path integrals asso- 



ciated with physically interesting quantities are difficult 
to evaluate because the action S in the exponent contains 
interaction terms (terms of third or higher order in the 
fields and their derivatives). In Casimir physics, on the 
other hand, the action is not more than quadratic in 0, 



and the difficulty in evaluating expressions like ( 39 1 stems 



instead from the challenge of implementing the implicit 
constraint on the functional integration measure, arising 
from the boundary conditions and indicated by the [• • • ]c 
notation in ( 39 ) . 

The innovation of Bordag |48] and of Li and Kar- 
dar |49j was to represent these constraints explicitly 
through the use of functional S functions. If the boundary 
conditions on (j) may be expressed as the vanishing of a set 



of quantities {La(j)}, where {La} will generally be some 
family of linear integrodifferential operators indexed by a 
discrete or continuous label a, then the constrained path 
integral may be written in the form 



m.) 



I?0„(x)l e-^[^"^«"l 
ic 

I?(/.„(x)n5(L,0)e- 



(42) 



where now the functional integration over 0„ is uncon- 
strained. A particularly convenient representation for the 
one-dimensional Dirac 5 function is 



Siu) 



dX 

2^' 



AXu 



(43) 



where we may think of A as a Lagrange multiplier enforc- 
ing the constraint that u vanish. Inserting one copy of 
( 43 1 for each 5 function in the product in ( 42 1 yields 



dXa 



ziU = JvM^)Jll^ 



^S[0,.:6>]+»E„A»io 



The final step is to evaluate the unconstrained integral 
over 0; since the exponent is quadratic in 0, this can 
be done exactly using standard techniques of Gaussian 
integration, yielding an expression of the form 



2{^n 



W/n 



dXr, 



-S^'ffAo} 



(44) 



[where {#} is a constant that cancels in the ratios in (40 



41 1]. The constrained functional integral over the field 
is thus replaced by a new integral over the set of Lagrange 
multipliers {Aq}, with an effective action S'^^ describing 
interactions mediated by the original fluctuating field (p. 



Representation of Boundary Conditions 



Equation MM makes clear that the practical conve- 
nience of path-integral Casimir computations is entirely 
determined by the choice of the Lagrange multipliers 
{Aq} and the complexity of their effective action S"^; 
these, in turn, depend on the details of the boundary 
conditions imposed on the fluctuating field. For a given 
physical situation there may be multiple ways to express 
the boundary conditions, each of which will generally 
lead to a distinct expression for the integral in (44 1. Ul- 



timately, of course, all choices must lead to equivalent 
results, but different choices may exhibit significant dif- 
ferences in computational complexity and in the range of 
geometries that can be efficiently treated. Several differ- 
ent representations of boundary conditions and Lagrange 
multipliers have appeared in the literature to date. 

The original work of Bordag et al. [48] considered QED 
in the presence of superconducting boundaries, with the 
boundary conditions taken to be the vanishing of the 
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normal components of the dual field-strength tensor; in 
the notation of the previous section, Lx0 = n'^FJlui^): 
and the set of Lagrange multipliers {Ax} constitutes a 
three-component auxiliary field defined on the bounding 
surfaces. The method is applicable to the computation 
of electromagnetic Casimir energies, but the treatment 
was restricted to the case of parallel planar boundaries. 

Li and Kardar [351 \EU\ considered a scalar field sat- 
isfying Dirichlet or Neumann boundary conditions on a 
prescribed boundary manifold. Here again the boundary 
conditions amount to the vanishing of a local operator 
applied to 0, Lx0 = 0(x) (Dirichlet) or ix0 — \d(l>/dn\x 
(Neumann), and we have one Lagrange multiplier A(x) 
for each point on the boundary manifold. In this case it is 
tempting to interpret A(x) as a scalar source density, con- 
fined to the boundary surfaces and with a self-interaction 
induced by the fluctuations of the cj) field. This formu- 
lation was capable, in principle, of handling arbitrarily- 
shaped boundary surfaces, but was restricted to the case 
of scalar fields. 

The technique of Refs. [49l [50] was subsequently re- 
formulated [Sni HH [53] in a way that allowed extension 
to the case of the electromagnetic field. Whereas the 
original formulation imposed a local form of the bound- 
ary conditions — and took the Lagrange multipliers A(x) 
to be local surface quantities — the revised formulation 
abandons the surface-source picture in favor of an al- 
ternative viewpoint emphasizing incoming and outgoing 
electromagnetic waves. This approach associates one La- 
grange multiplier A^ to each multipole term in a multi- 
pole expansion of the EM field, with the choice of mul- 
tipole basis (spherical, cyhndrical, etc.) governed by the 
symmetries of the problem; the effective action 5'°^ then 
describes the interactions among multipoles. 

The virtue of multipole expansions is that, for certain 
geometries, a small number of multipole coefficients may 
suffice to solve many problems of interest to high accu- 
racy. This has long been understood in domains such as 
electrostatics and scattering theory, and in recent years 
has been impressively demonstrated in the Casimir con- 
text as well 36, 52, 53J, where multipole expansions have 
been used to obtain rapidly convergent and even ana- 
lytically tractable series for Casimir energies in certain 
special geometries. The trick, of course, is that the very 
definition of the multipoles already encodes a significant 
amount of information about the geometry, thus requir- 
ing relatively little additional work to pin down what 
more remains to be said in any particular situation. 

But this blessing becomes a curse when we seek a uni- 
fied formalism capable of treating all geometries on an 
equal footing. The very geometric specificity of the mul- 
tipole description, which so streamlines the treatment of 
compatible or nearly-compatible geometries, has the op- 
posite effect of complicating the treatment of incompat- 
ible geometries; thus, whereas a basis of spherical mul- 
tipoles might allow highly efficient treatment of inter- 
acting spheres or nearly-spherical bodies, it would be a 
particularly unwieldy choice for the description of cylin- 



ders, tetrahedra, or parallelepipeds. Of course, for each 
new geometric configuration we could simply redefine our 
multipole expansion and correspondingly re-implement 
the full arsenal of computational machinery (a strategy 
pursued for a dizzying array of geometries in Ref. |36|). 
but such a procedure contradicts the spirit of a single, 
general-purpose scheme into which we simply plug an ar- 
bitrary experimental geometry and turn a crank. 

Instead, the goal of designing a more general-purpose 
implementation of the path-integral Casimir paradigm 
leads us to seek a representation of the boundary con- 
ditions that, while inevitably less efficient than spherical 
multipoles for spheres (or cylindrical multipoles for cylin- 
ders, or ...) has the flexibility to handle all manner of sur- 
faces within a single computational framework. This is 
one motivation for the fluctuating-surface-current (FSC) 
approach to Casimir computations, whose path-integral 
derivation we now discuss. 



B. Path-Integral Derivation of the FSC Casimir 
Formulas 



As noted above, key features of the path-integral treat- 
ment presented here include an unusual choice of action 
for the electromagnetic field and the introduction of sur- 
face currents as Lagrange multipliers constraining the 
photon field. After discussing these points in Sections 



IV B 1 and IV B 2 respectively, we show in Section IV B 3 
how together they allow us to evaluate the constrained 
path integral for the Casimir energy to obtain equation 



1. Euclidean Lagrangian for the electromagnetic field 

The usual (Minkowski-space) Lagrangian for the elec- 
tromagnetic field is 



S^J^ Jd^Ciiu,^), 



/:(w,x) = ^(6(c^,x)|E(c^,x)p-^(c^,x)|H(L.,x^'2 



Rewriting E and H in terms of the four-vector potential 
A^, integrating by parts, and rotating to Euclidean space 
via the prescription {w, A", ^"*} -> {i^, iA", iA°*} yields 
a Euclidean action density of the form 

Ce{^, x) =£i%iEl ( _ f A- A« - zCA'>*d.A^ 



+ 



1 



i^A'*d,A°-^A^*d,d,A° 



A'*djdjA' ~ A"d,djA^ 



or, introducing a convenient matrix-vector notation. 
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(45) 



(46) 



^dA 






The new four- vector field A^ defined by ( 46 ) will be the 
field over which we path-integrate, and equation (45) 



is almost, but not quite, the quantity that enters into 
the exponent of the constrained path-integral expression 
( 39 ). To complete the story, we must add a Fadeev-Popov 
gauge fixing term, which we do in analogy to the usual 
QED procedure [54] by simply displacing the coefficient 
of 5?2 term in (45) away from unity to ensure that the 
matrix in square brackets has no zero eigenvalues. Our 
final Euclidean action is 



Boundary conditions enforced by surface- current 
Lagrange multipliers 



In a scattering geometry consisting of one or more ho- 
mogeneous bodies bodies embedded in a homogeneous 
medium, the boundary conditions on the electromagnetic 
field are simply that the tangential E and H fields be con- 
tinuous across all material boundaries: if x is a point on 
the surface of an object, then we require 



/:E(C,x) = ^^^s)i(e)-(i 



2>2(0 



A"" 



^U 



(47) 



where the Faddeev-Popov gauge-choice parameter app 
may be chosen to have any finite value and is absent from 
all final physical predictions, as will be explicitly verified 
below [see equations ( 64p5 )]. Following the general pro- 
cedure reviewed in Section |IV A[ we can now write the 
Casimir energy at inverse temperature /3 in the form 



t''- 


E'"(x) - 


-E°"*(x) 


= 


(50a) 


b"- 


H-(x) - 


-H°"'(x) 


= 


(50b) 



where {E, H}'"'™* are the fields evaluated just inside 
and just outside the object surface at x, and where t" 
(a e {1, 2}) are vectors tangent to the surface at x [Figure 
l3](a)]. In terms of the modified four- vector potential A^^ 
these conditions may be written in the form 



_, 00 



Z{fi) 



2oo(/3,en) 



(48) 



i^{L-:°(x)-L^r(x)}^^(x)=0 (51a) 

i^{L--(x)-L-r(x)}^'^(x)=0 (51b) 



ZW,0 



VA" 



-|/A-S)(C)--4,dx 



(49) 



with the notation [• • • ]c indicating that the functional in- 
tegration ranges only over field configurations that satisfy 
the boundary conditions in the presence of our interact- 
ing material objects. 



where L'^^'' and L"'*" are differential operators that oper- 
ate on A'^ to yield the components of the E and H fields 
in region r. (We are here using a shorthand in which the 
A^ fields in the different regions, A'^"^ and y^^'°"*, are 
abbreviated simply as A^ and pulled outside the braces.) 
In a homogeneous region with spatially constant relative 
permittivity and permeability e(^,x) = e'"(^), /i(^, x) = 
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(a) 



{E,H} 



out 




out A out 




(c) 



FIG. 3. Enforcing boundary conditions via surface-current Lagrange multipliers, (a) Consider a single point x on the surface of 
an object in a Casimir geometry. The boundary conditions at x, which must be satisfied in the constrained path integral ( |49[ l, 
are that the tangential components of the E and H fields be continuous as we pass from inside to outside the object [equation 
(50l]; here {t^,t^} are vectors tangent to the surface at x. (b) We rewrite the boundary conditions in terms of differential 
operators L^'^'' operating on the A field, and we introduce Lagrange multipliers {K^ ,K^ , N^ , N^} to enforce the boundary 
conditions at x; specifically, K^ , K^ enforce the tangential E-field continuity at x, while 7V^,7V^ enforce the tangential H-field 
continuity [equation (|53[)]. (c) Repeating this procedure for all points on the object surface, we obtain Lagrange multiplier 



fields K(x),N(x), which have an obvious interpretation as the electric and magnetic surface currents of Figure [I] Integrating 
the photon field out of the path integral then yields an effective action describing the interactions of these surface currents 
[equations (581, (71 1, and (72|], leading ultimately to our fluctuating-surface-current formulas for the Casimir energy. 



/i''(^), the L operators take the form 



/ 



LE,r ^ 



1 






\ V^' 



d. li 





o\ 


dy 


^i 





d, 





.i) 



(52a) 



1 



/O -d, dy \ 



a^ -d^ 
\o -dy d^ / 



(52b) 



Equations (51 1 are a set of four boundary conditions for 



each point x on the surfaces of the material bodies in our 



geometry; in the language of Section IV A these are our 
constraints La4>, and to each constraint we now associate 
a Lagrange multiplier. We use the symbols K°'{x.) and 
iV°(x) (a = 1,2), respectively, to denote the Lagrange 
multipliers associated with constraints (51a,) and (51 d) 
at the single point x [Figure [31(b)] . Then the S functions 
that enforce the boundary conditions (51 1 at x are 



Ej"(x)-Ep(x) 



Hi|'(x)-Hr(x) 



dK 



{2ny 



>^ iK^-[Lj^.."-LE--']^^'(x) 



(53a) 



(27r)2 

(53b) 



where we may think of {Kx, Nx} = ELiI^''' ^"1*° as 
vectors in the tangent space to the boundary surface at x. 
Aggregating the corresponding d functions for all points 
on the surface of a single object, we obtain functional 
(5-functions, 



VK{x)e 



^/aoK(x)-[LE 



1.A-(x)dx (54^) 



I?N(x)e 



'/aoN(x)-[L^ 



'].A'"(x)dx 



(54b) 



where the integral in the exponent is over the surface 
do of an object in our geometry, and where the func- 
tional integrations / VK, J VN extend over all possible 
tangential vector fields on dO. 

Since K and N are tangential vector fields on dO that 
enforce the continuity of the tangential electric and mag- 
netic fields, respectively, it is tempting to interpret these 
quantities as electric and magnetic surface current densi- 
ties, and with their introduction our path-integral formal- 
ism begins to exhibit the first glimmers of resemblance to 
the surface-integral-equation picture reviewed in Section 

El 



3. Evaluation of the Constrained Path Integral 

In general we will have one copy of the functional 6- 
functions ( 54 ) for the surface of each object in our geome- 
try. Let {Kr,Nr} denote the Lagrange-multiplier distri- 
butions on the surface of the rth object; the constrained 
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path integral then reads 



VA'' 



■ J ATiAdx 



./nz>K,PN,/p^-|.-l/— ,55, 



with r,„ denoting integration over the surface of object 
r, and with the path-integration over A in the second 
hne now unconstrained [compare equation (42)]. This is 
just a standard Gaussian functional integral, which we 



proceed to evaluate using standard techniques [531 [5S] . 

To this end, it is convenient to think of breaking up the 
functional integration over A^ into separate integrations 
over the fields in each object and in the exterior region, 



fvA''= fvAi^YlVAf, 



where A!^ is the field in the interior of region r and A'^ is 
the field in the exterior region. The matrix S?(x), which 
depends on x through e and /i, is constant in each region 
due to the piecewise homogeneity of the geometry, while 
the operators L"^"''' only operate on the fields in region 
r. The functional integral becomes 



-f /^^ ^,.3,.^. -lE./v, Ar-Sr-A. 



2(/3,0 - fjlVKrVNr JVA'^, \{-DA^^[t 

r r 



with Jy denoting volume integration over the interior of 
region r. Now performing the Gaussian functional inte- 
grations over the fields Al^ immediately yields an expres- 
sion of the form ( 44 ) : 



2(/3,C) = {#}yn^K,I?N,e- 



1 ccff 



(57) 



where {#} is an unimportant constant that cancels upon 
taking the ratio in ( 48 1 (and which will be omitted from 
the equations below), and where the effective action for 
the surface currents. 



S° 



n 



N, 



^e 



{K.,N.}ti 



(58) 



contains terms describing both the self-interactions and 
the mutual interactions of currents on the object surfaces, 
as we now discuss. 

Consider first what happens when we integrate A!^ out 
of ([56]). Because the exponent of (56) couples A!^ only 



to Kr,Nr (and not to currents on other objects r' ^ r), 
integrating out Al;i yields an effective action involving 
only Kr, N^: 






The effective action Sr describes the self-interactions of 
electric and magnetic currents on dOr mediated by vir- 
tual photons propagating through the interior of object 
r. More precisely, we have 



K^,N^ (59) 

^ f dx f dx'l K^(x) • 7'=^'(x,x') • K^(x') 

+K,(x).7™-(x,x')-N,(x') 

+N,(x).7"--(x,x')-K,(x') 

+N,(x)-7™-(x,x')-N,(x' 



with the components of the tensor kernels given by (P,Q G 
{E,H}) 



P,)"^-! TQ,r 



PQ.r _ T'^^'^Tl^ T^' 



(60) 



We will see presently that the 7^'^^'' matrices here turn 
out to be nothing but the usual dyadic Green's tensors 
ppQ.r foj- tiie homogeneous medium inside object O^. 

To see this, it is easiest to perform the matrix multipli- 
cations of equation (60 1 in momentum space, where we 
have 
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and [cf. equation! 47 1 



S>7' = 



2>,M - 1 - 



a 



2) 



r,2 



IkP 



M 



app 



|k|2 



S'r,2 
S)r.2 



(61) 
(62) 



where 1 is the 4x4 unit matrix, {e'',/i''} are the (spatiaUy constant) permittivity and permeabihty of object r at 
imaginary frequency (,, k = \JfT [i"' ■ £,, and the momentum-space form of the £)r.2 matrix is 



2)r.2 = -- 



/ ^2 


Kkx 


Kky 


nkz \ 


Kkx 


kl 


iVj" ri"u 


kxkz 
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'^y^x 
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rZytZ^ 


\ nkz 


kzkx 




^i 1 



(The passage from the second to the third line of equation (62 1 is a standard algebraic manipulation in quantum field 
theory; see, e.g., equation (9.58) of Ref. [51]. The fact that we obtain such a concise form for the inverse of the matrix 
in square brackets in (61 1 is due to the fact that 2)^,2 is a rank-one matrix, and is known in numerical analysis as the 



Sherman- Morris on formula 



•) 



We now note the crucial fact that the second term in (62) makes no contribution to the effective action for the 



surface currents. Indeed, the contribution of this term to the 7 kernels in (60) involves triple matrix products of the 
form L^'*" • 'Z>r,2 ■ (L''''")"'"; but an explicit calculation reveals that (c = l/y/T^fF) 
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and the other three possible L • 3? • L products also vanish identically: 

L^'"^ • 2)^,2 • (L^'^'-)T = L"'''' • 2)^,2 • (L^'"")^ = L"''- • 2)^,2 • (L"^'^)^ 



0. 



(63) 



(64) 



(65) 



This furnishes the promised demonstration that the gauge-choice parameter ctpp makes no appearance in the effective 
action (59), thus explicitly confirming the gauge invariance of our procedure. 



Having verified that only the first term in (60) contributes to the 7 kernels in (60), these kernels are now easy to 
evaluate. First, 



-yEE.r ^ 
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But a quick comparison with the momentum-space forms of the dyadic Green's functions in Appendix \K\ reveals this 
to be nothing but ^ times the electric-electric dyadic Green's function for Or, i-e. 



or, transforming back to real space, 



7""''^(e;k) = e-r-^'''(C;k) 



7""''^(?;x,x') = C-r-^''-(C;x,x') 



and the connection of our formalism to SIE methodology begins to come into even sharper relief. 

Next, 
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E.r -i /T M,rNT 

—ickx i^ 

—icky «^ 

—ickz i^ 
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-k, 
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(66) 

(67) 

(68) 



thy Kx 



and again comparing with Appendix [K\ reveals that we have simply 



7 



^•r'= 



Having established the obvious pattern, it is now a short step to confirm that the remaining two cases of the 7 kernel 
in ( 60 ) are simply ^57, 



7" 
7" 






ME,r 
MM.r 



(69) 
(70) 



and equation (59), the portion of the effective action for surface currents that arises from integrating the photon field 

K,(x) 



in the interior of object Or out of (|56|), reads 



K,,N, 



'" L'"'VN.(x) 



"pEE.r "pEM,r 
"pA-iE.r "pMM,r 



K,(x') 
N,(x') 



(71) 



Next, we consider integrating the photon field in the exterior region (^J?) out of equation ([56|. Although the compu- 
tations proceed exactly as before, the the resulting contribution to the effective action is slightly more complicated. 
Because the exponent of (56 1 contains terms that couple A'^ to the currents on all object surfaces (unlike A'^, which 



couples only to surface currents on the single object Or), the result of integrating out A'^ will be an effective ac- 
tion describing the interactions of surface currents on all object surfaces mediated by exchange of virtual photons 
propagating through the external medium: 



{K,,N,} 



2^ 



dx / rfx' , . 



K,(x) 



"pEE,e "pEM,e 
"pME,e "pMM,e 



K,,(x') 
N,,(x') 



(72) 



Combining (71 1 and (72), the full path integral (57 1 now reads, for the particular case of two objects, 

T 



Z(/3,0 = y PK.PN.exp J -^y 





(73) 



with 



"pEE.l I "pEE.e 
■pME.l _|_ "pME.e 
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■pME,e "pMM,2 _|_ "pMM,e 
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Note that the quantity T ■ (j^) in ( 73 1 is nothing but the left-hand side of equation ( 14 ) . 

Having elucidated the structure or the effective action for surface currents, the remainder of our derivation is now 
straightforward. To evaluate the functional integral over sources in (73 1, we approximate K and N as expansions in a 



finite set of tangential basis functions representing surface currents flowing on the surface of the interacting Casimir 
objects, just as we did in equation (13): 



K(x)-^fc„f«(x), N(x) = -^^ 



We now insert these expansions into ( |73[ ) and approximate the infinite-dimensional integrals over K and N as finite- 
dimensional integrals over the k^, and Hq, coefhcients: 



X>K,I?N, 



J jYl dka n '^'^o 



where J7, the Jacobian of the variable transformation, is an unimportant constant that cancels upon taking the ratio 
in ( 48 ) (and which will not be written out in the equations below) . Equation ( 73 ) now becomes simply 



2(/3,e)=/n*«n^-"^"*^"^''^*'^'^"^' 



a finite-dimensional Gaussian integral which we evaluate immediately to obtain 



{#}-[detM(0 



-1/2 



(74) 



where M(^) is nothing but the SIE matrix discussed in Section ffn (and, once again, |# | is just an irrelevant constant 
into which the ^/2/? prefactor in the exponent disappears). Now finally inserting ( |74[ ) into (491 and (48) leads 
immediately back to our FSC formulae ([l]), and our derivation is complete. 



V. EQUALITY OF THE PARTIAL TRACES 

An important practical simplification of the FSC for- 
mulae follows from the structure of the BEM matrices. 
Recall from (II]) that the quantity that enters into the 
FSC expression for the Casimir force is 



Tr <^ M" 



dM 






dM, 



fSa 



dvi 



(75) 



(with a similar expression for the torque). We will show 
that the trace in ( 75 ) neatly splits into two equal sub- 
traces, and thus that to compute the full trace we need 
only sum a subset of the diagonal elements of the matrix 
in curly brackets (and double the result). In practice this 
reduces the computational expense of the trace compu- 
tation by a factor of two or greater. 

To understand the physical intuition behind this sim- 
plification, note that the sum in ( 75 ) runs over the ele- 



ments of our basis of surface-current expansion functions, 
which includes functions defined on the surfaces of each 
of the objects in our Casimir geometry. In evaluating 
the portion of this sum contributed by basis functions 
defined on a single object surface, we are in effect com- 
puting a sort of surface integral over the surface of that 
object. Intuitively we might expect that, to compute the 
force on one object, it would suffice to evaluate this sur- 
face integral over the surface of that object alone, or over 



r 



the surfaces of all other objects, but that we need not do 
both. As we now show, this physical expectation is born 
out by the mathematics; to compute the full trace in ( 75 ) 



we need only sum the contributions of basis functions on 
the surface of the object on which we are computing the 
force — and double the result. 

As before, let the objects in our Casimir geometry be 
labeled Or, r — 1,2,- ■■ , with Oi the object on which we 
are computing the Casimir force, and let Dr be the di- 
mension of the subblock of the M matrix corresponding 
to object Or (that is, Dr is the number of surface-current 
basis functions defined on the surface of 0^)- The dimen- 
sion of the full matrix is D = J2 Dr ■ 

A rigid displacement of object Oi leaves unchanged 
the interactions between all pairs of basis functions save 
those pairs in which precisely one basis function lives on 
Oi. We can thus split the sum in (75) into two pieces: 



E^.-; 



dM, 



aP 



dr^ 
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^E E ^a 

a=lp=Di+l 
D Di 

■ E E^c.^ 

Q=Di-|-l 0=1 




(76) 



The first piece on the RHS here is the sum of the first 
Di diagonal elements of the matrix M"^ • ^^, while the 
second piece is the sum of the remaining D—Di elements. 
But from the fact that M is a symmetric matrix (Map = 
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Mi3a) it now follows that the two pieces here are equal, 
and thus to compute the full trace we need only sum 
the first Di diagonal elements of M^^ • ^^ (or the latter 
D — Di elements, if they are fewer) and double the result, 



TrM" 



dM 
dvi 



2E E ^-J 



a=ip=Di- 



dJVL 



a/3 



dvi 



In practice, a convenient way to evaluate this quantity 
is to LU-factorize the matrix M, solve the linear sys- 
tems M • Xm — Bm where the vectors H-m are the first 
Di columns of ^^, then extract and sum the mth el- 
ements of the vectors B^ and double the result. The 
equality of the partial traces then ensures that this op- 
eration requires just Di linear solves, in contrast to the 
full D solves that would be required in the absence of the 
simplification. 



VI. APPLICATIONS 



A. Casimir Forces between Metallic Spheres and 
Cubes 



To validate our new method and demonstrate its flex- 
ibility, we first calculate the Casimir force between pairs 
of metallic particles of spherical and cubical shapes (Fig- 
ure |4|. 

For the sphere-sphere case, the Casimir force may 
be computed using scattering-matrix methods based on 
a spherical-wave decomposition of the electromagnetic 
field |52l I58j , with the Casimir energy expressed in terms 
of interactions among waves labeled by the usual spher- 
ical indices {£,m}. At large values of the sphere-sphere 
separation, only waves with small values oi £ are relevant, 
and in this regime Ref. [52] obtained an asymptotic power 
series for the force, whose first four terms we have plot- 
ted in Figure I4] (lower dashed curve). At smaller values 
of the sphere-sphere separation, the sum over spherical 
waves may be evaluated numerically, as is done in Figure 
[4] (hollow blue circles) with waves up to ^ = 40 retained 
for each of the two possible polarizations, corresponding 
to 2{£ + 1)^ = 3362 basis functions for each sphere. 

To perform the calculation using FSC techniques, we 
discretize the surfaces of the objects into small surface 
patches and expand surface currents using the localized 
basis functions of Figure [flA] in this case, the number of 
basis functions retained in the description of the surface 
currents is 2976 for each sphere. As illustrated in Figure 
|4J the FSC calculation (solid red circles) reproduces the 
results of the spherical-wave calculation. Note that our 
choice of surface-mesh basis functions allows us to con- 
centrate more degrees of freedom in the regions where 
we expect the surface source densities to be most rapidly 
varying — namely, the regions of each sphere that most 
closely approach the other sphere — while simultaneously 
using a relatively coarse-grained representation of weakly 



interacting regions. Thus already for the simple sphere- 
sphere geometry our method begins to exhibit practical 
advantages over spherical-wave-basis methods, in which 
the resolution can only be increased globally rather than 
locally. 

The FSC method really comes into its own in treating 
objects that cannot be efficiently described by analyti- 
cal Maxwell solutions, such as the case of two cubes. To 
handle such a geometry using scattering-matrix methods, 
we would be forced either to expand fields in and around 
the cubes in spherical waves — an approach which would 
require retaining basis functions up to inordinately large 
values of £ except in the long-distance limit [51] — or to 
reformulate the method in a basis of Maxwell solutions 
for cubical scatterers, for which analytical expressions 
are not available. In contrast, the FSC approach han- 
dles the geometry with no more effort than is required 
for the sphere-sphere case (inset and solid blue circles in 
Figure |4]). 

An additional limitation of the spherical-wave ap- 
proach is that it is inherently restricted to separation 
distances large enough that the interacting objects may 
be enclosed in non-touching, nonoverlapping spheres. 
For the cube-cube geometry pictured in Figure Q, 
this would render the method inapplicable for distances 
L/R < 1.47, excluding much of the range plotted in the 
figure. 



B. Repulsion of an anisotropic nanoparticle from a 
square aperture in a thin metallic plate 

Reference [5] considered the Casimir force on an elon- 
gated nanoparticle above an aperture in a thin metallic 
plate and predicted a regime in which the force on the 
nanoparticle is repulsive. This work, as well as subse- 
quent investigations of similar phenomena [501 WT\ , con- 
sidered only circular apertures; here we investigate the 
case of a square aperture (Figure p]) . We consider the 
Casimir force on a cylindrical nanoparticle (lower inset), 
whose axis coincides with the axis of a thin (20 nm) plate 
with a 1 /xm square aperture, as a function of the dis- 
tance Z between the vertical center of the nanoparti- 
cle and the center of the thin plate. Both cylinder and 
plate are made of real (lossy) gold, described by a rel- 
ative dielectric function e(^) = 1-1- Wp/[^(C + 7)] with 
{wp, 7} = {1.37- 10l^ 5.23- 1013} rad/s. When the center 
of the nanoparticle is vertically aligned with the center of 
the plate {Z = 0), the z-directed Casimir force vanishes 
by symmetry; as the nanoparticle is displaced slightly 
in the positive z direction it experiences first a repulsive 
Casimir force (shaded region of plot) which peaks near 
the value of Z at which the nanoparticle first lies entirely 
above the plate. As the nanoparticle is raised further 
above the plate, the repulsive force decreases in magni- 
tude and eventually crosses over into an attractive force 
whose magnitude decays at large Z (unshared region of 
the plot). 
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FIG. 4. Casimir force between PEC spheres and between PEC cubes. The filled green circles indicate sphere-sphere data 
computed using the FSC method described in this paper, while the hollow blue circles indicate sphere-sphere data computed 
using a numerical implementation of the scattering-matrix method of Ref. ^5] • (The lower dashed-line curve indicates the first 
four terms of the asymptotic series for the sphere-sphere Casimir force reported in Ref. [52].) The filled red squares indicate 
cube-cube data computed using the FSC method described in this paper; for this geometry, scattering-matrix methods and 
indeed almost all existing Casimir methods would be unwieldy or impossible to apply. The upper dashed-line curve indicates 
the proximity-force approximation (PFA) to the cube-cube force, -FppA = 4-R^ ■ ttjtjt^ for cubes of face area 4R^. 



The fact that the square-hole geometry reproduces 
the repulsion phenomenon observed in the circular-hole 
case is not particularly surprising, but we note that this 
geometry exhibits several features which would make 
it prohibitively expensive if not outright impossible to 
treat using any other Casimir method of which we are 
aware. In particular, the interpenetrating nature of the 
nanoparticle-plate configuration would immediately ren- 
der most scattering-matrix methods inapplicable. More- 
over, the drastic distance in length scales between the 
tiny (20 nni diameter) cylinder and the relatively large 
(5 fj,ia radius) plate, together with the absence of rota- 
tional symmetry, would pose severe challenges to finite- 



difference methods. In contrast, the FSC method eas- 
ily accomodates nonuniform surface meshes, allowing 
us to describe the nanoparticle surface as a union of 
appropriately-sized surface patches even as we use much 
larger patches for the plate surface. 



VII. CONCLUSIONS 

The fluctuating-surface-current (FSC) approach to 
Casimir physics, together with other recently-developed 
Casimir methods such as the finite-difference ap- 
proaches [3ni [SB [S21 [S3] J constitutes an advance in the 
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FIG. 5. Casimir force on an elongated cylindrical nanoparticle above a square aperture in a thin plate. The upper inset (not 
to scale) shows the particle-plate geometry, while the lower inset is a close-up view of the surface mesh used to represent the 
particle. The axis of the nanoparticle (the z-axis) is perpendicular to the plane of the plate; the center of the nanoparticle lies 
a distance Z above the center of the plate. In the regime < Z < 330 nm, the z-directed force on the nanoparticle is positive, 
i.e. the particle is repelled from the plate; the shaded portion of the graph indicates this repulsive-force regime. For Z > 185 
the nanoparticle lies entirely above the plate and the force is thus unambiguously repulsive. Both the nanoparticle and the 
plate are made of finite-conductivity gold (see text). 



development of Casimir algorithms that mirrors an ear- 
lier evolution in computational electromagnetism (EM). 
In the latter field, traditional special-function approaches 
such as Mie's method for spherical scattering began to be 
complemented in the 1960s and 70s by a host of numer- 
ical techniques that expanded the range of geometries 
that could be accurately and efficiently treated. Today, 
numerical techniques for general geometries coexist with 
geometry-specific special-function approaches to consti- 
tute a rich and varied arsenal of computational EM tech- 
niques appropriate for almost any conceivable situation. 

In the future, we expect a similar situation to pre- 
vail in the Casimir field. General-purpose methods such 



as the one presented here will not replace scattering- 
matrix Casimir methods, any more than finite-difference 
Maxwell solvers have replaced the theory of Mie scat- 
tering; instead, methods such as our FSC technique 
will coexist with scattering-matrix methods, augmenting 
the toolbox of available methods available for predicting 
Casimir interactions across the gamut of geometric and 
material configurations. 

Among general-geometry numerical Casimir methods, 
the FSC approach is unique in obtaining compact deter- 
minant and trace formulas for Casimir quantities. All 
other general-purpose numerical methods rely on numer- 
ical evaluation of a surface integral for the Casimir force 
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and torque (and on even more unwieldy numerical volume 
integrations for the Casimir energy). For this reason, we 
expect the FSC approach to be the most efficient nu- 
merical Casimir method for the piecewise-homogeneous 
material configurations to which it applies. 

In addition to the practical usefulness of the FSC 
Casimir formulas, the two independent derivations that 
we have provided in this paper contain a number of 
theoretical innovations that we expect to find broader 
application. In particular, in the stress-tensor deriva- 
tion we stated and proved a new integral identity in- 
volving the homogeneous dyadic Green's functions of 
Maxwell's equations (Appendix IbI), while in the path- 
integral derivation we introduced a new type of Lagrange 
multiplier to constrain the functional integration over the 
photon field. We hope the latter technique will prove 
to be a generally useful tool in quantum field theory; 
one possible application beyond the Casimir realm is the 
boundary dependence of entanglement entropy in the 
electromagnetic field, a subject recently addressed for 
scalar fields f64]. 

What challenges lie ahead for FSC Casimir computa- 
tions? If the evolution of computational Casimir physics 
continues to mimic that of computational electromag- 
netism, an obvious next step will be the development of 
fast solvers [55] — algorithms that exploit physical insight 
to reduce the computational complexity of matrix ma- 
nipulations in equations like M from O(iV^) to a more 
tractable scaling such as O(A^logA^), where N, the di- 
mension of the matrix, is the number of surface-current 
expansion functions {fa} retained in equations Q and 
(13). Although a number of algorithms are known for 



evaluating matrix-vector products involving SIE matri- 
ces in 0{N\ogN) time |5S], the question of how best to 
exploit such algorithms to evaluate the determinant and 
trace in equations (fTl) is nontrivial. To date, all FSC 
Casimir calculations have considered problems of moder- 
ate complexity {N < 10'*), for which dense-direct linear- 
algebra solvers are adequate; however, future problems 
may require going beyond this regime, in which case fast 
solvers will be essential. 

A separate challenge is to apply FSC techniques to the 
calculation of Casimir forces out of thermal equilibrium, 
as well as to the closely related problem of near-field ra- 
diative heat transfer. These problems, which represent 
a logical next step in the study of fiuctuation-induced 
forces beyond the equilibrium Casimir case, have be- 
gun very recently to be investigated using scattering- 
matrix [261 1661 167[ and numerical 68J methods. A hy- 
brid technique combining scattering-matrix ideas with 
SIE-based numerical calculations was proposed in 8J , but 
was restricted to a cylindrical-wave basis; to date there 
has been no basis-independent surface-current formula- 
tion of nonequilibrium fluctuation problems. Could the 
FSC technique presented in this paper be modified to ap- 
ply to these problems? If so, what physical insight would 
the surface-current formulation lend, and what improve- 
ments could be achieved in the efficiency of practical cal- 



culations? 

In this paper we have said almost nothing about the 
practical challenges inherent in concrete numerical imple- 
mentations of the FSC formulas. Among these are (a) 
how best to choose the surface-current basis functions 
{fa} for a given geometry, (b) how to evaluate the mul- 
tidimensional integrals that enter into the elements of the 
M matrix [such as equations ([7|]; (c) how to compute 
the matrix determinant, inverse, and trace in equations 
(fTl), and (d) how to evaluate the imaginary- frequency 
integrations in equations (fTl), as well as the Matsubara 
sums in their finite-temperature analogues. All of these 
challenges will be addressed in subsequent publications. 
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Appendix A: Homogeneous and Scattering Dyadic 
Green's Functions 



For reference, and to define our notation, we col- 
lect here some well-known results concerning the dyadic 
Green's functions of classical electromagnetic theory [131 

fSSl- 

In the presence of known volume densities of electric 
and magnetic current J(x),M(x) at a fixed imaginary 
frequency ^, the components of the electric and magnetic 
fields are given by linear convolution relations of the form 

i?,:(x) = y'{r^-(^;x,x')J,(x') +rf;,"(e;x,x')Af,(x')}dx' 

ff.(x) ^ y'{r^.-(e;x,xV,(x') + r™(e;x,x')A/,(x')}dx'. 

These relations define the four dyadic Green's functions 

r. 

In an infinite homogeneous medium with spatially con- 
stant relative permeability and permittivity e(^;x) = 
e''(^), /i(^;x) = fJ.^{£,), the four T functions may be ex- 
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pressed in terms of just two tensors: 

r™''^(C, r, r') = -Z^Z^'kJ' G(K^ r - r') 
r"'=''^(C,r,r') = «:''C(K'',r-r') 
r™''^($, r, r') = -k" C(K^ r - r') 



r™''^(C,r,r') 



Z^Z 



:G(«:^r-r') 



(Ala) 
(Alb) 
(Ale) 

(Aid) 



^0 



Z'' 



k'' = VA'oM'^eoe'' • C 



where G, sometimes referred to as the "photon Green's 
function," is the solution to the equation 



V X V X +K^ 

and C is defined by 



G(K;r) =(S(r)l 



C = -V X G. 

K 



(A2) 



(A3) 



(Note that G and C have dimensions of inverse length, 
while the F functions have dimensions of field/surface 
current density; for example, F""^ has dimensions of mag- 
netic field / electric surface current density.) 

Explicit expressions for the components of G and C 
are 



Gij - 



S,j - ^d,dj 



Go, 



rv 



StjkdkGo (A4) 



where Gq is the scalar Green's function for the Helmholtz 
equation. 



Go(«:;r-r') = 



^-K|r-r'| 



(A5) 



47r|r - r'l 
which satisfies 

V^ - K^j Go{k; r - r') - S{r - r'). 

Momentum-Space Representations 

The momentum-space decomposition of the scalar 



Green's function (A5) is 



with 



Go{^;r)^ I ^^Go(«;k)e^''"- 



Go(K;k) 



(In what follows we will generally omit the ^ designa- 
tion, relying on context to differentiate between real- and 
momentum-space functions.) 



The momentum-space version of (A4) reads 
S. 



rZ-i Ki j 



Gij — 



or, in matrix format. 



Go 



Gij = — SijkkkGo (A6) 

til 



G(k) = 



1 

K2(K2 + |k|2 



K^ 

^2 
k2 



(A7a) 



C(k) = 



K(K2 + |k|2) 






(A7b) 



Ky Kx 



Scattering Dyadic Green's Functions 

In a general inhomogeneous region, the dyadic Green's 
functions may be expressed as the sum of two terms, 

F--(e; X, x') = F---(^; x - x') + ^""(C; x, x') (A8) 

(and similarly for the other three F functions); here F"^"^" 
is the homogeneous DGF for an infinite medium with 
constant e^ , ff set equal to their values at x, and Q^^ is 
the scattering part of the DGF, which describes the fields 
scattered from the inhomogeneities in the geometry. The 
first term in ( A8 1 is singular as x — > x' , but the sec- 



ond term is perfectly well-defined in that limit and is the 
quantity that enters into the fluctuation-dissipation ex- 
pressions for the spectral density of fluctuations in prod- 
ucts of field components, as discussed in Section |lTT| 



Appendix B: Proof of Integral Identities 

In this Appendix we state and prove a new integral 
identity that underlies the stress-tensor derivation of the 
FSC formulae presented in Section [llTl To our knowledge, 
this integral identity is new and is stated, proved, and 
used for the first time in this paper, although an iden- 
tity bearing at least a superficial resemblance appears in 
equation A. 6 of Ref. 70j. 

In SectionlTTIj we introduced a three- index integral ker- 
nel X defined by 

I,fe,(r,r') = K' /v:4(r,r',x)n^(x)dx (Bl) 

where the integration is over a closed surface in space 
(a topological two-sphere) C and the integrand contains 
products of factors of dyadic Green's functions: 



Va{v,v',^) 






5iA 



Gl y-^2 y-^1 y-^2 j^ 



'Bk^lB' 



2 '^Bk'-'IB- 

(B2) 
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FIG. 6. Schematic summary of the integral identity (B5l. proved in Appendix B We consider a closed bounding surface C 



and choose two points r and r', which may lie bot h inside (left panel), both outside (center panel), or on opposite sides (right 
panel) of C. We write an expression [equation (B4|] that involves products of dyadic Green's functions (DGFs), one connecting 
r to a point x on C and a second connecting x to r'. Then we evaluate the surface integral of this expression as x ranges over 
all of C. The answer we obtain depends on the relative positioning of r and r' with respect to C: If the two points lie both 
inside or both outside C, the surface integral vanishes, while if the two points lie on opposite sides of C then the surface integral 
yields the derivative of a DGF connecting r to r'. 



Here the G and C dyadics are those defined by equations 
( |A2[ ) and (A3); we use capital Roman letters (A,i3, • • • ) 
to denote contracted indices, and we are using a short- 
hand notation in which k arguments are suppressed and 



spatial arguments are replaced by superscripts. 



Gj^- = G'ij(K,x-r), 



Gi 



Gy(K,r'-x). 



We also defined a symmetrized version of I ; 

Iifc;(r,r') =Zife;(r,r') +Zi;fe(r',r) 



= n^ 



FA(x)nA(x)dx 



with 



^^(x) = g^g: 









(B3) 



(B4) 



The goal of this Appendix is to demonstrate that, by appealing to the defining properties of the G and C dyadics, 
the surface integral in (B3) can he evaluated in closed form, with the result (depicted schematically in Figure pi) 



Tiki(Y,r') 



0, 



dr. 



if both r, r' lie inside C 



■Gkiir — r') if r lies inside and r' lies outside C 



d_ 
dvi 



(B5) 



[0, 



Gfe/(r — r') if r lies outside and r' lies inside C 
if both r, r' lie outside C. 



[As far as we can tell, the symmetrization in equation B3 is necessary to achieve the compact form of equation 



( |B5[ ); our attempts to evaluate the non-symmetrized X in concise form were unsuccessful.] 

For the purposes of this Appendix it will be convenient to work in length units such that n — \. With this convention, 
the G and C dyadics are related to the scalar Green's function for the Helmholtz equation according to 



G, 



did J Go, 



Oi. 



-£ijkdkGo, 



(B6) 
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and Go satisfies 



OaBa - 1 



Go(r) = 5{v). 



(B7) 



The evaluation of the surface integral ( B3 1 now proceeds in several stages 



1. Apply Divergence Theorem 



The first step is to recast the surface integral in (B3| as a volume integral over the volume bounded by C, 



VAnAdA^ I dAVAdV 
c Jv 



WdV 



C = dV 



where we put W = OaVa- 

Each of the six terms in V a contains two factors and hence contributes two terms to W (by the chain rule for 
differentiation). Terms of the form dACiA vanish; to each of the remaining (nonvanishing) terms we assign a label, as 
tabulated in Table II Note that derivatives of Gfj enter with a minus sign, because we are differentiating with respect 
to X, which enters the argument of Gf; with a minus sign. 

TABLE II. Terms in W = OaVa- 



Label 


Term 


W^ 


dAGik GiA 


W^ 


- Glk dAGfA 


W^ 


^aG^a Gil 


w^ 


- GIa OaGu 


w^ 


— diGkA GiA 



Label 


Term 


W"^ 


GkA (^iGlA 


w-' 


- dACfk CIa 


W^ 


[cIa] [dAGf^ 


w^ 


OiCAk C[A 


w^" 


^ GAk diCiA 



2. Treatment of W^ and W^ 



We first consider the terms labeled W^ and W^ in Table pl| Starting with the first of these, we have 



W^ 



G,i 



ik 



dAGfA 



Expand the second factor using ( B6 ) 



Gik 



di{l~dAdA)Gl 



Apply (B7) 



= + Gj, diSi2) 
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Integrate by parts: 



diGik 



m 



(B8) 



where ~ means "equivalent as long as we are underneath the volume-integration sign," and ^(2) is shorthand for 
5{r' — x). (Note that the minus sign coming from the integration by parts is cancelled by the minus sign coming from 
the fact that x enters 2 with a minus sign, as noted above.) 
By analogous operations, we find 



W"" 



5{1) 



dkGii 



(B9) 



where (5(1) is shorthand for (5(x — r). We now set aside results (B8) and (B9) for future use. 



3. Treatment of Remaining Terms 

a. Rewrite in terms of Go 



Turning next to the remaining 8 terms in Table 1, we begin by using (B6) to rewrite everything in terms of the 
scalar Green's function: 



W^ + W^ + W^ + w^ 



^aG^}. — diGf. 



kA 



G 



lA 



GkA 



^i^iA ~ (^aG^, 



SikdAGl - SkAdiGl 



SiaG^ - OiOaG^ 



+ 



SkAGQ — dkdAGQ SiaOiGq — SuOaGq 



W 



W^ + W^ + W^° 



■^^ikB^lAC 
SAkB^lAC 



OaObGI 

^^^BGl 



dcG^o 
dcGl 



' ^kAB^ilC 
^AkBSlAC 



ObGq 
ObGI 



dAdcGl 
didcGl 



We can trade Levi-Civita symbols for Kronecker deltas using the standard identity 
sabc^def 



= 5al 



SreS, 



BEOCF — OBFOCE 



SbfSc 



Sai 



SrfS, 



BFOCD — OBDOCF 



SroSc 



Sai 



SrdS, 



BDOCE — OBEOCD 



SreSc 



We then find 



w'^ + W^ + W^ + w^° 



OaOaGI 


[SudkG 


2 




SkAG 


2 




- 6 


ikdlGg - 


SkAGl 


OaOaG^ 


-Su 


dkdAGl 


OaG^ 


+ S.,k 


OaGI 


OiOaG^ 




- 


d. 


Gl] 


'dkdi 


Gi 


+ 


dkdiG 


1 




d^G 


2 

n • 







(BIO) 



(Bll) 



6. Label Individual Terms 



To proceed, we now assign a new label to each separate term in (BIO) and (Bll ): 



16 



w^ + w^ + w^ + w^ + w'^ + W^ + W^ + w^° = ^ ^" 



n=l 



TABLE III. Terms in X. 
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Label 


Term 


X^ 


—Ski diGo Go 


X^ 


[9,Gj] [dkdiG^'^ 


x^ 


&ik diGo Go 


x^ 


-S,k[dAGl'^\^dAdiG'i] 


x^ 


—5ii Go dkGo 


x« 


Ski Go diGo 


x^ 


5ii[dAdkGl'\[dAG^'\ 


x« 


- [dkdiGl\ [acg] 



Label 


Term 


x"" 


-Ski [dAdAGl] [a.Gg] 


^10 


+9»! [a^^AGj] [dkG^'] 


X" 


+Ski [ftGj] [^A^AGg] 


j^l2 


-<5»fc [diGl^ [OaOaG^] 


^13 


—Sii dAdkGo BaGo 


j^l4 


+S,k [BaGI^ [didAGi\ 


^15 


- [diGi\ [dkdiGl] 


^16 


[dkdiGl^ [^^Gf^ 



c. Recombine Terms 



We first note the obvious cancellations: 

X'^ + X^^ = X^ + X" = 
X'^ + X^^ = X^ + X^^ = 

Summing and appropriately recombining the remaining terms, we find 

^1 + ^3 + X^^ + X^^ = [SkAGl - S,kdiGl\ [dAdAG^ - Gg 

and 



X^+X^ + X^+X^° 



dAdAGQ — Gq 



5iidkG^ — SkidiGQ 



From ( B6 1 it follows that 



SkidiGo — SikdiGo — diGki — diGik 



Using this and (B7), we can rewrite the previous two equations in the form 

X^ + X^ + X" + X^^ = \d,Gli - diGJ^ \6{2) 

and 



(B12) 



X^ + X"" + X"" + X^° = (5(1) dkGl-d,Gli 



(B13) 
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4. Final Steps 



Finally, we combine equations (|B8|), (|B9|), (|B12|) and (|B13|) to obtain 

1-ikl = 



I iw'^ + W^ +X^ + X^ + X^^ + X^'^ + X^ +X^ + X^ ^ X^°\ dy 



Writing out the function arguments, this reads 



Iiki{r,r') 



diGki{x - r) (5(r' - x) - d{x - r) diGki{r' - x 



dx. 



(B14) 



We now proceed on a case-by-case basis depending on the positions of r, r'. 

• First, if r and r' both lie outside the bounding surface C, then neither 6 function contributes and we have 1 = 0. 

• If r lies inside C while r' lies outside C, then only the second S function contributes, and we find 

lM{r,r')^~d,Gki{r'-r) 
= +d,Gkiir~r') 

• If r lies outside C while r' lies inside C, then only the first S function contributes, and we find 

Im{r,r')^d,Gki{r'-r) 

^-d,Gki{r-r'). 

• Finally, if both r and r' lie inside C then both 6 functions contribute, their contributions cancel, and we find 

i = o! 



The result (B5) is thus established, and our proof is complete. 



The J Kernel 



Th e de velop ment of Section |IIIC| also makes referen ce to a version of the I kernel defined in analogy to equations 
(iBll) and (|B2|, but with the ''GG - GG" structure of (IbI replaced by a ''CG + GG" structure: 



Jifci(r,r') = K^ * YA(r,r',x)nA(x)dx 

^A(r, r ,x) — G^j^Gij^ Y^^k^lB + '^ik^lA Y^Bk^iB- 



(B15) 
(B16) 



We also define a symmetrized version defined in analogy to (B3): 



J^kl{r,r') = J^kl{r,Y') + Jak{v' ,v). 

In Appendix [A] we noted that the curl operation takes G into C and C into — G. [Technically, in the latter case there 
is an additional 5 function, which we neglect for reasons discussed in conjunction with equations (70) above.] With 
this observation we see that (B16) is obtained from (B2) simply by taking the curl with respect to the k index, and 



thus that the J kernel is the result of the same operation applied to the I kernel: 

1 



Jikii'T.v') = -efe^B9Al,B;(r,r') 



= < 



0, 



9r, 



if both r, r' lie inside C 



Cfei(r — r') if r lies inside and r' lies outside C 



d_ 



0, 



Cfe;(r — r') if r lies outside and r' lies inside C 
if both r, r' lie outside C. 
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